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Spatially explicit and temporally continuous estimates of photosynthesis will be of great importance for in- 
creasing our understanding of and ultimately closing the terrestrial carbon cycle. Current capabilities to 
model photosynthesis, however, are limited by accurate enough representations of the complexity of the un- 
derlying biochemical processes and the numerous environmental constraints imposed upon plant primary 
production. A potentially powerful alternative to model photosynthesis through these indirect observations 
is the use of multi-angular satellite data to infer light-use efficiency (e) directly from spectral reflectance 
properties in connection with canopy shadow fractions. Hall et al. (this issue) introduced a new approach 
for predicting gross ecosystem production that would allow the use of such observations in a data assimila- 
tion mode to obtain spatially explicit variations in 8 from infrequent polar-orbiting satellite observations, 
while meteorological data are used to account for the more dynamic responses of 8 to variations in environmen- 
tal conditions caused by changes in weather and illumination. In this second part of the study we implement and 
validate the approach of Hall et al. (this issue) across an ecologically diverse array of eight flux-tower sites in 
North America using data acquired from the Compact High Resolution Imaging Spectroradiometer (CHRIS) 
and eddy-flux observations. Our results show significantly enhanced estimates of s and therefore cumulative 
gross ecosystem production (GEP) over the course of one year at all examined sites. We also demonstrate that 
s is greatly heterogeneous even across small study areas. Data assimilation and direct inference of GEP from 
space using a new, proposed sensor could therefore be a significant step towards closing the terrestrial carbon 
cycle. 

© 2012 Elsevier Inc. All rights reserved. 


1. Introduction 

Global models of plant photosynthesis are limited in their ability 
to infer the spatial and temporal heterogeneity of photosynthetic 
light-use efficiency (s), the efficiency with which plants can use 
absorbed radiation energy to produce biomass. This e is driven by nu- 
merous environmental constraints affecting the photochemical reac- 
tion chain, and varies greatly in space and time. Our limited 
understanding of the underlying biochemical processes (Field & 
Mooney, 1986) and difficulties in obtaining the driving variables 
globally are major limitations to current approaches modeling photo- 
synthesis (Turner et al., 2003). 


* Corresponding author. Tel.: + 1 301 286 8597: fax: + 1 301 614 6695. 
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A potentially powerful alternative to deriving e from environmen- 
tal constraints is to directly infer the status of the xanthophyll cycle, a 
biochemical mechanism that balances light use and absorption in 
higher plants (Demmig-Adams & Adams, 1996). Under conditions 
where light use efficiency is limited by factors other than light, the 
xanthophyll cycle pigment violaxanthin is rapidly converted via in- 
termediate antheraxanthin to zeaxanthin, both of which have bio- 
chemical structures that allow the dissipation of excessive radiation 
energy as heat. This pigment conversion can be quantified by the pho- 
tochemical reflectance index (PRI), a narrow waveband index that 
uses a xanthophyll specific absorption band at 531 nm, however, con- 
founding background and bi-directional reflectance effects have ham- 
pered its use for almost two decades (Coops et al., 2010). Recent 
progress using multi-angular observations(Hall et al., 2008), has 
allowed us, for the first time to infer s across a large range of forested 
ecosystems from space using one functional relationship, the first 
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derivative of PRI with respect to shadow fractions (ct s ), hereafter re- 
ferred to as PRI'. Theoretical (Hall et al., 2011, 2008) and empirical 
work (Hilker et al„ 2011, 2010a) has given this multi-angular method 
a solid foundation as generic, physically based measure of light-use 
efficiency: PRI' is insensitive to background reflectance and reflec- 
tance of non-photosynthetically active material, because for single 
leaf scattering, a valid assumption for wavelengths in the visible 
bands, Hall et al. (2008) showed that PRI cannot change its value 
with shadow fraction unless the value of one of its bands changes as 
a physiological response to the degree of leaf illumination. Hilker et 
al. (2011) implemented these theoretical considerations over a 
range of different study sites using multi-angle satellite data. The re- 
sults are summarized in Fig. 1 and illustrate the capacity of this meth- 
od to infer e across a large range of at least temperate and boreal 
ecosystems from space. As a result, this new technique could bypass 
much of the difficulties experienced in modeling e from environmen- 
tal constraints. 

In the first part of this study (Hall et al., this issue) we described a 
simple approach for modeling e from data assimilation of PRI' and 
meteorological constraints. In the forward mode, our model down 
regulates a spatially explicit maximum convergence efficiency (E opt ) 
using a non-linear, multivariate response function of photosyntheti- 
cally active radiation (PAR), temperature (T) and water vapor pres- 
sure deficit (D) to obtain temporally continuous estimates of GEP at 
a 30 m spatial resolution. This multivariate function accounts for the 
variability in e due to the highly dynamic changes in the environment 
in between satellite overpasses. However, since e can also be ob- 
served directly during a satellite overpass, these spaceborne observa- 
tions can be used together with corresponding measures of PAR, T 
and D during the spacecrafts overpass time to invert the model and 
infer a pixel specific e opt . This e opt can be re-determined every few 
satellite overpasses, and as a result, can be used to describe more 
inert changes in the conditions affecting e (such as soil nutrient status 
and edaphic water stresses for instance), thereby reducing much of 
the complexity of modeling GEP(Hall et al., this issue). 

In this second paper we evaluate our new approach across a num- 
ber of temperate forest sites using observations from the Compact 
High Resolution Imaging Spectroradiometer (CHRIS), a demonstra- 
tion sensor of the European Space Agency (ESA) aboard the Platform 
for On Board Autonomy (Proba). First, meteorological response 
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Fig. 1. Relationship between PRI' as observed from CHRIS/Proba imagery and EC mea- 
sured £ for DF-49, Harvard Forest, HJP1975, HJP2002, Howland Forest, Northern Old 
Black Spruce, Old Jack Pine and Southern Old Aspen. All observations have been collect- 
ed between 2001 and 2009. The study demonstrated that even over a large different 
canopies, a consistent relationship was can be established between PRI' and e, hence 
PRI' serves as a generic tool to measure e remotely. 

From Hilker et al. (2011) Journal of Geophysical Research, Vol. 116, G03014, Reproduced 
by permission of American Geophysical Union. 


functions are evaluated using eddy flux observations made at eight 
different tower sites throughout North America. Our T and D func- 
tions are compared to response functions of the MODIS GEP algo- 
rithm (MOD17) (Heinsch et al., 2006; Running et al., 2004; Zhao et 
al., 2005), a commonly used technique to assessing global carbon up- 
take from earth observing satellites. Also, we use CHRIS/Proba data to 
investigate the spatial heterogeneity of s opt across and within the dif- 
ferent study sites used in this paper. 


2. Methods 

2.1. Study sites 

Eight different sites were selected to evaluate our model. The se- 
lected sites span a broad range of temperate forest types despite a 
limited range due to the availability of simultaneous eddy-flux 
tower observations and CHRIS/Proba satellite data. Table 1 provides 
a description of the study locations and data availability. The struc- 
ture of the different sites ranged from recently harvested, regenerat- 
ing Jack pine ( Pinus banksiana Lamb.) dominated forest (HJP2002) to 
mature, Douglas-Fir ( Pseudotsuga menziesii var menziesii (Mirb.)) 
dominated (e.g. DF-49) and Aspen (Populus tremuloides Michx) dom- 
inated (e.g. SOA) stands (see Fig. 2 for in situ photographs of all sites). 
For all but one site (SOA), CHRIS/Proba data were available in con- 
junction with eddy flux observations and have been collected be- 
tween 2001 and 2009. 


2.2. Eddy flux and meteorological data 

Canopy GEP was determined from eddy-covariance (EC) measure- 
ments acquired from the data archive of the Canadian Carbon Pro- 
gram (CCP) for the Canadian sites and Ameriflux for the two US 
sites (Table 1). Net ecosystem exchange (NEE) was determined as 
the sum of the half-hourly fluxes of C0 2 and the rate of change in 
C0 2 storage in the air column between ground and EC- 
measurement level (Barr et al., 2004), fluxes for the Harvard site 
were available as hourly observations (Horii et al., 2004). 

Incident and reflected PAR [pmol m~ 2 s _1 ] was measured from 
upward and downward looking quantum sensors above and below 
the canopy and the fraction of PAR absorbed by the canopy (f PAR ) 
was derived at each site from the incident and reflected total PAR 
measured above and below the canopy (pi(d) and p 2 (0J), the effec- 
tive leaf area index (L e ), and the solar zenith angle (0) at the time of 
measurement (Chen, 1996). GEP was determined as the difference 
between NEE and daytime ecosystem respiration (R D ) (Humphreys 
et al., 2006) as provided by the respective flux network. R D was calcu- 
lated using the annual exponential relationship between nighttime 
NEE and soil temperature at 5-cm depth after applying a logarithmic 
transformation to correct for heteroscedasticity (Black et al., 1996). 
Light-use efficiency was derived from (Monteith, 1972, 1977) 


GPP 

PAR x f p AR 


( 1 ) 


More detailed descriptions on processing of the eddy-covariance 
data can be found in (Humphreys et al„ 2006), a comprehensive re- 
view on Fluxnet procedures and processing of EC-data is provided 
in (Baldocchi, 2003). 

Above-canopy air temperature (T) and relative humidity were de- 
rived from temperature and humidity probes housed in aspirated 
shields (Humphreys et al., 2006). Atmospheric water vapor pressure 
deficit (D) was computed from T and relative humidity (Buck, 1981) 
and PAR was derived directly from the incoming PAR sensors. 
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Table 1 

Site description and available CHRIS/Proba scenes. The overpass time varies by up to several minutes per site. All times are given in GMT. 
Site description is adapted from Hilker et al., 2011. 


Site, Reference 

Lat (°)/ 
Long (°) 

Elev 

(m) 

Dominant species 

Leaf area 
index 

Age 

(years) 

Height 

(m) 

Annual 
mean temp. 

("O 

# of avail. 

CHRIS- scenes /approx, 
overpass time 

DF-49 (Morgenstern et al., 2004) 

-125.334 

49.867 

340 

Pseudotsuga menziesii, Thuja 
plicata, Tsuga heterophylla 

7.1 

60 

35 

8.1 

6 

18:45 

Harvard (Munger et al., 1996; Staebler & 
Fitzjarrald, 2005) 

-72.171 

42.537 

340 

Quercus rubra, Acer rubrum, Betula 
lenta, Pinus strobus, Tsuga canad 

3.4 

80 

23 

8.3 

2 

16:01 

HJP1975 (Amiro et al., 2006; Chen et al., 

2006; Schwalm et al., 2006; Zha et al., 2009) 

-104.645 

53.876 

570 

Pinus banksiana 

1.4 

35 

6 

0.4 

9 

18:09 

HJP2002 (Amiro et al., 2006; Chen et al., 2006; 
Kidston et al., 2010; Schwalm et al., 2006) 

-104.649 

53.908 

560 

Pinus banksiana 

0.9 

8 

0.1 

0.4 

9 

18:09 

Howland (Hollinger et al., 2004; 
Xiao et al., 2005) 

-68.740 

45.204 

60 

Picea rubens, Tsuga canadensis 

5.3 

109 

20 

6.7 

3 

15:18 

NOBS (Bergeron et al., 2007) 

-98.481 

55.880 

259 

Picea mariana 

4.8 

160 

9 

-4.4 

6 

18:14 

OJP (Amiro et al., 2006; Chen et al., 2006; 
Schwalm et al., 2006) 

-104.692 

53.916 

579 

Pinus banksiana 

2 

91 

13 

0.4 

9 

18:09 

SOA (Barr et al., 2004) 

106.198 
- 52.629 

600 

Populus tremuloides 
Corylus comuta 

2.1 

83 

22 

0.4 

0 


2.3. Shadow fractions 

Estimation of canopy shading (a s ) is critical for accurate modeling 
of canopy light-use efficiency, as sunlit parts of the canopy are more 
likely to be exposed to excessive radiation energy than shaded 


vegetation elements (Hilker et al., 2008b). For instance, Hall et al. 
(2008) showed that canopy level measurements of the photochemi- 
cal reflectance index (PRI) are strongly dependent on a s , and that 
the directional changes observed in PRI at a given half hour interval 
can be attributed almost entirely to changes in ot s (Hall et al„ 2008). 



Fig. 2. In situ photographs and overview map of the 8 research sites presented in this study. The sites are DF-49 (A), Harvard Forest (B), HJP1 975(C), HJP2002 (D), Howland Forest 
(E), Northern Old Black Spruce (F), Old Jack Pine (G) and Southern Old Aspen (H). 
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In previous studies (Hall et al., 2008; Hilker et al„ 2008a, 2010a) we 
have used airborne light detection and ranging (LiDAR) to estimate 
canopy shading. While this approach yields a good approximation 
of mutual shading of individual stands, LiDAR data were not available 
for all of the sites. As a result, we determined a s at all flux-tower sites 
using a geometric optical model (Li & Strahler, 1985), which esti- 
mates mutual canopy shading using measures of tree and crown di- 
mensions and estimates tree density. Those measurements were 
available from previous Held measurements (Hilker et al., 2010b) 
and published values (Amiro et al„ 2006; Bergeron et al., 2007; 
Chen et al., 2006; Schwalm et al., 2006; Staebler & Fitzjarrald, 2005; 
Xiao et al., 2005) for the different flux tower sites. 

Canopy shading effectively increases the proportion of diffuse sky 
radiation by blocking direct sunlight from parts of the canopy. The 
magnitude of this effect depends on the fraction of diffuse sky radia- 
tion and is largest under direct illumination conditions, while non- 
existent under completely cloudy skies. While not all of the sites 
used in these studies provided direct measurements of the diffuse 
sky radiation, the proportion of diffuse radiation can be approximated 
when comparing the measured incident PAR to modeled clear sky 
conditions as a function of the solar zenith angle (Hilker et al., 
2009a). Clear sky solar irradiance was modeled from maximum PAR 
observations obtained for a range of zenith angles at each site and 
sky conditions were estimated as the fraction of actually measured 
irradiance to the maximum solar potential for that zenith angle at a 
given point in time. For details see (Hilker et al., 2009a). 

2.4. CHRIS/Proba imagery 

CHRIS is an imaging spectrometer with a 615 km sun- 
synchronous orbit and an orbital repeat cycle of approximately 
7 days. The CHRIS/Proba configuration permits along-track narrow- 
band spectrometric observations of PR1 of up to five angles ( + 55°, 
+ 36°, 0", —36°, —55°) at a maximum spatial resolution between 
18 m or 34 m depending on the data acquisition mode. This data is ac- 
quired nearly simultaneously within each overpass during which 
stand level e may be considered constant. CHRIS/Proba images ac- 
quired between 2001 and 2009 were obtained for 8 of the test sites 
(no data were available for the SOA site) from ESA’s online archive 
(https://oa-es.eo.esa.int/ra/). The image size of CHRIS/Proba scenes 
is approximately 25 x 25 km at nadir. 

CHRIS data collected in modes 1 and 3 were used, as they provided 
the appropriate PRI wavebands at around 531 and 570 nm. The satel- 
lite observations were pre-processed using the V1SAT tool of the 
European Space Agency (Gomez-Chova et al„ 2008), to convert satel- 
lite measured radiance to top of atmosphere reflectance and screen 
the images for clouds (Thuillier et al., 2003). A two step geo- 
rectification algorithm was applied (Ma et al., 2010) and CHRIS/ 
Proba satellite images were registered with respect to Landsat obser- 
vations of the same locations (Hilker et al., 2011). First, common 
ground control points (GCPs) between Landsat and CHRIS/Proba 
were automatically identified using a scale-invariant feature trans- 
form (SIFT) (Lowe, 2004). The network of these initial GCPs was 
then enhanced in a second step using a normalized cross correlation 
(NCC) approach. 

Satellite observed reflectance depends on two main parameters, 
aerosol optical thickness (AOT) and surface reflectance (SR). To cor- 
rect for the impact of atmospheric effects, commonly used, pixel 
based algorithms of atmospheric correction assume a Lambertian sur- 
face reflectance. Hilker et al. (2009a, 2009b) showed that this as- 
sumption directly contradicts the multi-angle effects observed in 
PRI and consequently, no meaningful multi-angular observations of 
PRI can be obtained from single orbit atmospheric corrections 
(Hilker et al., 2009b). This challenge could potentially be overcome 
when using a time series approach for which the retrieval of aerosol 
optical thickness does not require these simplifying assumptions 


(Hilker et al., 2009b; Lyapustin et al., 2007); however, no such algo- 
rithm currently exists for CHRIS/Proba. Using radiative transfer theo- 
ry, it can be shown that atmospheric path scattering will cause PRI to 
increase with as hence depress PRI' and the estimate of [epsilon] 
(Hall et al., 2011). PRF will be affected by the effects of different 
path lengths through the atmosphere, but not the absolute change 
in PRI due to lack of atmospheric correction. The limitations due the 
lack of atmospheric correction are acknowledged, the calculated 
error (Hall et al., 2011) should, however, still allow the demonstra- 
tion of the proposed model from top of atmosphere reflectance. 

PRI was computed from CHRIS/Proba imagery as the normalized 
difference of CHRIS bands 4 (529 nm, Bandwidth: 12.9 nm) and 6 
(569 nm, Bandwidth: 14.1 nm) for images acquired in CHRIS Mode 
3 (all sites except for the southern BOREAS region), and band 11 
(532 nm, Bandwidth 13.4 nm) and 15 (573 nm, Bandwidth 9.6 nm) 
for images acquired in CHRIS Mode 1 (Hilker et al., 2011 ). PRI and cor- 
responding a s were computed for each pixel of a multi-angular image 
stack acquired during one overpass (Hilker et al., 2011). Please note 
that the slope of PRI with respect to canopy shading (a s ) is non linear, 
as demonstrated in Hall et al. (2011). However, this non-linearity is 
most important near the dark spot. For the purpose of this study, 
we have used a linear regression of PRI vs. a s to simplify the proces- 
sing, as the CHRIS/Proba data observes all scenes at roughly the same 
range of shadow fractions at all times, and these view angles, and 
Hilker et al. (201 1 ) showed that for this data the slope can be approx- 
imated as simple linear regression. 

2.5. Model inversion 

Inversion of our multivariate model (Hall et al., this issue) allows 
the parameterization of £ opt across the landscape from CHRIS/Proba 
reflectance. In the forward mode, the photosynthesis model intro- 
duced in this study provides estimates of light-use efficiency for a 
given set of landscape parameters, i.e. D, T, PAR and e opt . In the in- 
verse mode, the CHRIS/Proba inferred e can be used to derive e opt 
for every given pixel of a satellite scene, assuming that solar irradi- 
ance, D and T are constant across the scene (which is reasonable to 
a first approximation given the limited spatial extend of CHRIS/ 
Proba of approximately 20 x 20 km). 

Mathematically, model inversion is a non-linear minimization 
problem that can be solved through iterative adjustment of estimated 
a-priori inputs (Verstraete et al., 1996). Different optimization algo- 
rithms are available, based on cost functions to minimize the resid- 
uals between forward modeled and measured observations (here: 
canopy reflectance). In this study, we selected a trust-region- 
reflective algorithm based on the interior-reflective Newton method 
(Coleman & Li, 1996; Coleman et al., 2002). The search range of s opt 
was set to 0-4gCMJ _1 for all sites based on an extensive review 
by Schwalm et al. (2006). 

Technically, model inversion would also allow inference of the 
shape parameters introduced for the response functions of D, T and 
PAR. However, the limited availability of CHRIS/Proba data makes a 
dynamic inference of these parameters difficult. As a result, we deter- 
mined the parameters of the response functions empirically, by 
inverting the model using eddy-flux observations and assuming the 
shape parameters to be constant across all sites. 

3. Results 

The first section of the results presents a comparison between the 
response function introduced in Hall et al. (this issue) and EC-flux 
data. In order to evaluate the degree to which our multivariate, 
non-linear response functions improves the GEP estimates over the 
MODIS MODI 7 approach, we use the MODI 7 performance as a base- 
line. The second section shows the results from model inversion using 
the CHRIS/Proba data. 
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Fig. 3 shows the dependency of e on D and T and PAR across the eight 
different research sites selected for this study. The X-axis represents the 
temperature response curve for a given level of D, the Y-axis shows the 
D response curve for a given temperature level. The characteristic 
shapes of both response functions are clearly visible in the model (com- 
pare Hall et al., this issue; Figs. 3, 4 and 7). The shape of the three dimen- 
sional mesh-surface reflects the multi-variate nature of our model. For 
instance, when T is near the optimum, e tends to be higher for a given 
level of D, because one stress factor has less of an impact than if two 
stresses affect photosynthesis at the same time. The dots represent 


actual observations of T, D and corresponding e made by the eddy co- 
variance systems and averaged as hourly observations. The colors of 
the data points correspond to different levels of PAR (red = high, blue 
= low). As the physiology predicts, high levels of PAR yielded low levels 
of e. Again, down-regulation is more rapid, when one or both of the 
other limiting factors (T and D) are sub-optimal. For instance, higher 
PAR levels yielded high e values only when temperatures were near op- 
timum ( see for instance Figure B) ; otherwise e was low, even if PAR was 
low. The three dimensional mesh represents the fit of the physiological- 
ly based response functions to the data, the colors of the mesh-surface 
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Fig. 3. Dependency of s on D and T using hourly EC-flux data (only daily data were available in cases of B, C and G) for 8 different tower sites. The 3D surface is showing the bounding 
envelope, the dots are the actual observations. The color code of the mesh corresponds to the different levels of s, the color code of the observations corresponds to the level of PAR 
(dark = low, bright = high). 
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Fig. 4. Mean and standard deviation of minimum and maximum 8-day temperature, D and PAR at the 8 research sites and annual course of the T, D and PAR modifiers, shown in red. 
The T and D modifiers of the MODIS algorithm are shown in blue. The left column shows the temperature modifier, the central column represents the modifier for the vapor pres- 
sure deficit, and the right column shows the PAR modifier. 


highlight the different levels of 8. The response function envelop reflects 
the maximum s inferred for a given temperature and water vapor pres- 
sure deficit. Our model assumes that all further reduction in s is caused 
by increases in PAR. 

Fig. 4 shows a comparison between the meteorological response 
functions for D and T using the traditional approach and ours (Figs. 4 
and 5 in Hall et al., this issue). The response functions presented in 
Fig. 4 are based on biological rationale and their coefficients were es- 
timated by minimizing the difference between tower-measured GEP 
and model computed GEP. One year worth of data is shown for each 
site. The left column in Fig. 4 shows the minimum (gray) and 


maximum (black) temperature, here averaged over eight days for dis- 
play purposes. The error bars in all columns represent the standard de- 
viations around the means. The central column is showing variations 
in D (mean and standard deviation); the right column represents 8- 
day variations in PAR (mean and standard deviation). The correspond- 
ing response functions of the MODI 7 algorithm are shown in blue 
(right axis of the left and central columns), a response function of 
the model introduced in this study is shown in red. As it was shown 
in Fig. 6 of Hall et al. (this issue), 8 responses to T, D and PAR are not 
separable into multipliers because the 8 responses to T, D and PAR 
are interdependent. Consequently, the red lines in Fig. 4 show the 
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Fig. 5. Comparison between EC derived £ (daily averages) and daily averaged estimates of e determined using the MODI 7 approach for the 8 different flux tower sites. 


response of the e modifier to one of the three variables with the other 
two fixed, when all other stresses are at a minimum (smallest possible 
reductions). The response functions for the MOD17 model and the 
model introduced in this study were computed from hourly observa- 
tions of eddy flux data and then averaged over the eight day period 
for display purposes. 

A clear seasonal signal was found at each site, however, seasonal 
differences were particularly strong at the three BOREAS sites, whereas 
the climate variability was smallest at the Coastal Douglas-fir (DF-49) 
site. When using the MODI 7 algorithm, D limited photosynthesis 
only during short periods of the year, or not at all, as in case of the 
Harvard site. 

The response functions observed from the MODI 7 algorithm 
showed a largely temperature driven limitation to photosynthesis. 
During the beginning and the end of the season, the temperature 
responses of our model were similar to those shown for the MODI 7 
algorithm; however, our model yielded a significant reduction in e 
during mid summer, when temperatures were high. Our model also 
showed a much stronger reduction of the modeled photosynthesis 
due to D (driving stomatal conductance) when compared to the 
MODI 7 algorithm, due to its non-linear character. The PAR-based 


reduction in photosynthetic efficiency was comparatively higher at 
the deciduous sites and mixed than at the coniferous sites due to 
the large amount of shading observed at the coniferous sites, causing 
an effective reduction in direct incident solar radiation. The PAR 
based response was strongest at the HJP 2002 site, which experienced 
almost no shading as a recently harvested, regenerating site. 

Fig. 5 shows a comparison between tower-derived e and the effi- 
ciency modeled from environmental constraints using the MODI 7 al- 
gorithm across the different research sites. The e opt values were 
obtained from the MODI 7 lookup table for the respective biome 
type. While the model captured some of the seasonal patterns 
shown in the eddy-flux data, values saturated at E opt for notable 
parts of the growing season at most of the observed sites. The algo- 
rithm tended to overestimate e, especially during the onset of the 
growing season and later in fall. The model also overestimated the 
boreal coniferous sites (Figure D, G) for most of the year. The differ- 
ences in modeled s resulted in an overestimation of GEP for most 
sites, except DF-49. The RMSE between tower observations and mod- 
eled e values was 0.32 gCMJ' 1 (DF-49), 0.40 gCMJ' 1 (Harvard), 
0.38 gCMJ' 1 (Howland), 0.33 gCMJ' 1 (OJP), 0.38 gCMJ' 1 
(HJP1975), 0.32 gCMJ" 1 (HJP2002), 0.32 gCMJ' 1 (NOBS) and 
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Fig. 6. Comparison between EC measured hourly GEP (only daily data were available in cases of B, C and G) and hourly, daily and 8-day averaged estimates of GEP determined using 
the MOD1S approach for all flux tower sites. The y-axis shows the cumulative sum in GEP acquired for up to one year in length. 


0.25 gCMJ" 1 (SOA). Fig. 6 shows a comparison of cumulative GEP 
estimated over one year as derived from the flux tower observations 
and the MOD17 algorithm. In order to compare the effects of different 
averaging techniques, data is shown for estimates of e made from 
hourly, daily and 8-day observations of T Min and D. The results 
shown in Fig. 6 have been normalized by the number of observations 
made, in order to allow a direct comparison of the cumulative sums. 

The use of non-linear response functions in combination with dy- 
namic £ opt values inferred from satellite observations helped in in- 
creasing the accuracy of estimated e and GEP values considerably. 
Fig. 7 shows examples of £ opt values derived across the different 
study areas from model inversion of CHRIS/Proba observations. 
Three of the study sites (OJP, HJP 2002, and HJP1975) where all con- 
tained in one CHRIS/Proba scene (Fig. 7D). The size of the maps 
shown in Fig. 7 varies from site to site, as £ opt requires a stack of 
multi-angular CHRIS/Proba scenes to overlap in order to be able to 
infer PRF (see Hilker et al„ 2011 for details). The figure illustrates 


the spatial heterogeneity of the estimated £ opt values at each of the 
sites. No multiple inferences of £ opt in time were made at this point, 
as most of the sites were limited by the availability of sufficient satel- 
lite data, and as a result all satellite images during the seasons dis- 
played in Fig. 7 were used to derive £ opt . However, the potentials of 
inferring multiple E opt throughout the seasons are illustrated by the 
example of the southern BOREAS area and shown in Fig. 10 (This 
site had 7 useful satellite images available throughout 2006 and 
2007). Highest variability in e opt was found at the most structurally 
diverse landscapes, especially the heavily logged area around the 
DF-49 site, but also at the southern Boreas site. Only moderate differ- 
ences in CHRIS/Proba-inferred £ opt were found at the NOBS site, 
which is one of the more homogeneous sites under investigation as 
illustrated in Fig. 2. No CHRIS/Proba data were available for the SOA 
site. 

A comparison between tower-derived e and the efficiency mod- 
eled using the non-linear responses are shown in Fig. 8. The £ opt 
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Fig. 7. Examples of 8 opt values derived across the different study areas from model inversion of CHRIS/Proba observations. (A: DF-49, B: Harvard Forest, C: Howland Forest, D: OJP, 
HJP and HJP 2002, E: NOBS). 


values were obtained from CHRIS/Proba as presented in Fig. 7 and 
averaged over an area of 3 x 3 pixels around the tower to include 
about 80% of the daytime flux footprint (Hilker et al., 2008a). In 
case of the SOA site, the maximum observed e as determined by the 
eddy covariance system was used as a surrogate (top 1% of the 
data). The higher E opt values prevented the model from saturating 
thus allowing it to capture most of the variability in s throughout 
the season. The RMSE between tower observations and modeled e 
values was 0.24 g C MJ -1 (DF-49), 0.20 g C MJ _1 (Harvard), 0.23 g C 
Mp 1 (Howland), 0.17 gC Nip 1 ( OJP), 0.16 g C Mp 1 (HJP1975), 


0.12 g C MJ -1 (HJP2002), 0.15 gCMp 1 (NOBS) and 0.15 gCIVp 1 
(SOA). The non-linear response also increased the range in variability 
in e thereby allowing it to follow the EC observed values more closely. 
As a result, the new model yielded more realistic estimates of GEP 
across the 8 different research sites with the errors in cumulate GEP 
being notably reduced compared to Fig. 6 (Fig. 9). 

Fig. 10 shows a time series of £ opt values inferred for the example 
of the southern Boreas study area between 2006 and 2007. s opt was 
derived from model inversion using three consecutive CHRIS/Proba 
scenes in a “moving window” approach, the dates provided in 
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Fig. 8. Comparison between EC derived 8 (daily averages) and daily averaged estimates of e determined using the approach introduced in this study for the 8 flux tower sites. 


Fig. 10 are representing the date of the first of three consecutive im- 
ages. Please note that the spatial extend of the maps varied from 
time step to time step as the common area of overlap was different 
for each acquisition. The figure illustrates a consistent spatial pattern 
of E op t across the southern Boreas study area between mid 2006 and 
mid 2007. While the temporal frequency of observations is not 
enough to allow a comprehensive analysis of the seasonal effects on 
e opt , some seasonal patterns could be observed. For instance, the 
e opt value for the OJP tower decreased consistently from 
2.7gCMJ -1 in June 4, 2006 to 2.4gCMJ _1 in October 2006 and 
went back up to 2.7 g C MJ~ 1 in August 2007. Over the same time pe- 
riod, g opt values for the area around the HJP 2002 tower showed a 
slight, but consistent decrease from 2.1gCMJ~ 1 to l.OgCMJ^ 1 
and went back up to 2.3 g C MJ~ 1 in August 2007. 

4. Discussion 

This study introduced and verified a, physiologically grounded 
model to infer landscape-level photosynthesis from data assimilation 
of meteorological observations and multi-angular satellite images. 


The technique previously developed in Hall et al. (201 1 ) and Hilker 
et al. (2011) allowed us, for the first time, to observe instantaneous 
e and therefore GEP directly from space (Hilker et al., 2011) across a 
wide range of forested ecosystems. The intention of the present 
work was to demonstrate how such information could be used in a 
data assimilation scheme to improve global estimates of GEP. The 
lack of spatial heterogeneity in biome specific E opt values has previ- 
ously been identified as a significant drawback (Goulden et al., 
1996) of the MODIS MOD17 approach (Running et al., 2004). The 
use of multi-angular satellite observations may help to overcome 
some of these limitations, by 1 ) allowing a spatially explicit mapping 
of E opt across the landscape and 2) allowing E opt to vary over time, at 
least in the presence of sufficient spaceborne observations. 

The non-linear response functions based on Bernacchi et al. 
(2002) and Jarvis (1976), but coupled as multivariate functions, 
yielded a notable improvement over the previous techniques based 
on linear datasets, at least for the sites investigated in this study. 
The non-linear model showed better results in capturing the seasonal 
variability in s (Figs. 5 and 7), whereas the linear functions tended to 
overestimate e during the early and late season and saturated for 
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Fig. 9. Comparison between EC measured hourly GEP (only daily data were available in cases of B, C and G) and hourly, daily and 8-day averaged estimates of GEP determined using 
the approach introduced in this study for 8 different flux tower sites. The y-axis shows the cumulative sum in GEP acquired for up to one year in length. 


much of the summer season (Fig. 5) (compare also Turner et al., 
2003), as little or no limitations were imposed by either D or temper- 
ature. The use of the non-linear response functions yielded a sharper 
limitation of s optt as a result of D, T min and PAR. Additionally, our 
model response to above optimal temperature conditions resulted 
in further reductions during the summer across most of the sites 
(perhaps with the exception of the northern BOREAS study region, 
NOBS, see Fig. 4) where warmer temperatures may be advantageous. 
By down-regulating photosynthesis with PAR, our multivariate re- 
sponse function also had a notable impact on s, especially at the de- 
ciduous and clear-cut sites. The stronger response with respect to 
environmental stress factors allowed for a more realistic £ opt which 
was close to the maximum e observed by the eddy covariance data 
at the eight sites (compare the peaks of modeled and observed values 
in Fig. 3). While Fig. 4 gives an impression of the differences in the re- 
sponse functions used in both approaches, it should be noted that the 
responses introduced in this model are multivariate and 


consequently, the functions in Fig. 4 only represent the reductions 
which are imposed if all other stresses are zero (compare Fig. 3). A di- 
rect comparison between the reduction levels of both algorithms is 
therefore difficult and should be viewed with caution as no interac- 
tion between the photosynthesis limiting factors is allowed in the 
MOD1S approach. The differences in modeled, cumulative GEP as 
shown in Figs. 6 and 8 demonstrate the advantages of the non- 
linear, multivariate, physiologically-based responses. Specifically, 
the deviations in cumulative GEP are a direct result of the differences 
in which e, as the estimates of f par and PAR were identical for EC- 
measured, linearly and non-linearly obtained models. 

While the non-linear functions introduced in Hall et al. (this issue) 
yielded more realistic estimates of e and therefore GEP, they also 
added more complexity. For instance, the shape parameters of the 
curves were derived from eddy covariance observations using 
model inversion, which requires a certain level of availability of 
such data across the landscape. Such data are currently available 
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Fig. 10. Time series of s op t values derived from CHRIS/Proba scenes over the southern BOREAS study area. The dates (DOY) of observations were 2006/155 (A), 2006/197 (B), 2006/ 
240 (C), 2006/274 (D), 2006/300 (E), 2007/187 (F) and 2007/214 (G). 
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through flux tower networks, assuming that the functional form of 
the stress response function is identical for all vegetation types across 
a given area or biome (note however, the quantity is still driven by 
8 op t, which is different for every 30x30 m pixel). This assumption is 
reasonable to a first approximation, as these responses are mainly 
shaped by evolutionary processes and are therefore largely a result 
of climate adaptations. As a result, the functional form (not the mag- 
nitude) of the response functions can be assumed to be relatively in- 
variant, at least across climatically similar regions. Limitations are 
faced with respect to differences in elevation, as higher regions, 
even if in spatial proximity, can vary significantly in climate. Similar 
assumptions, however, are also made in case of the MOD1S MODI 7 al- 
gorithm as minimum and maximum T min and D levels are determined 
for each biome type (Heinsch et al., 2006; Running et al., 2004; Zhao 
et al., 2005). 

Arguably, one of the most significant advantages of the technique 
introduced in this study is the inference of a spatially explicit, tempo- 
rally dynamic e opt from the CHRIS/Proba imagery. The maps shown in 
Fig. 9 illustrate the heterogeneity in s opt , even across the relatively 
small study areas. These differences are largely a result of changes 
in canopy structure and illustrate the need for an enhanced mapping 
of photosynthesis. The heterogeneity was largest at the coniferous 
sites (especially the southern BOREAS, but also the DF-49 site) as e 
in conifers is heavily driven by shadow fractions (Hilker et al., 
2008b, 2010a). The high e opt values in Fig. 7a are due to a lack of var- 
iation in shadow fractions at clear cut areas around the DF-49 site and 
present a limitation of the approach introduced in Hilker et al. (201 1 ) 
and Hall et al. (2011): The satellite based technique is using the rela- 
tionship between PRI and a s to obtain a generic function of e. This re- 
quires a sensor to observe vegetation under a certain range of a s , 
which may not be available at less structured, herb and grass domi- 
nated areas. In such cases, it may be possible to base such observa- 
tions primarily on PRI as the structure of these ecotones is much 
simpler, however, further study will be required to address this 
issue and also to investigate its significance with respect to terrestrial 
carbon budgets at a landscape scale. 

The time series of e opt values which were derived from CHRIS/ 
Proba data acquired between June 2006 and July 2007 illustrates the 
potential for a temporally dynamic inference of s opt from satellite ob- 
servations. Little variation was found in e opt throughout the observed 
time period which could suggest that most of the seasonal variability 
at this site was driven by changes in D, T and PAR. However, these 
findings have to be interpreted with caution, as no off-season images 
were available and the temporal density of the existing time series is 
too sparse to show conclusive evidence. In addition, we required mul- 
tiple observation dates to infer s opt from model inversion (Hall et al., 
this issue; in this study we chose 3 consecutive dates). As a result, 
each of the maps shown in Fig. 10 represents 1 to 6 months of data 
thereby potentially smoothing some of the seasonal patterns that 
might exist. Due to data limitations, the observations used in this 
study, reflect relatively “normal" years, more variability, also in LAI, 
would be expected when observing extreme events, such as severe 
drought stresses and further research may be required to assess the 
algorithms performance under these circumstances. It is expected 
that e opt would be greatly reduced and further reductions due to D 
T and PAR would largely prevent photosynthesis. 

For model inversion, E opt was limited to 4 g C MJ' 1 based on pre- 
vious, tower obtained findings (Schwalm et al., 2006). While we 
found that this value was a good threshold for forested ecosystems 
around the tower sites, higher thresholds may be defined for vegeta- 
tion with the potential for higher s opt values. For instance, some of the 
harvested or disturbed areas in Fig. 7 A and D show e opt values of close 
to 4gCMJ -1 . Furthermore, D and T were assumed to be constant 
across a CHRIS/Proba scene. While this is probably reasonable to a 
first approximation, it should be noted that both variables will vary 
with vegetation type and also with elevation, and consequently, 


higher spatial resolution observations of temperature and humidity 
would be helpful to address this concern. 

One of the biggest limitations of the algorithm presented in this 
study is the lack of suitable satellite data. To date, CHRIS/Proba is 
the only satellite in orbit that acquires PRI data from multiple view 
angles along track. The instrument, does not obtain data globally 
and, as a result, images exist for a very limited number of pre- 
defined test sites. The images are also not freely available, which fur- 
ther restricts its use. This lack of observations may result in limited 
possibilities for inferring E opt at regular intervals throughout a grow- 
ing season. As a result, the model may under-represent the seasonal 
effects of photosynthesis (as shown around DOY 200 in Fig. 8C and 
G). Hilker et al. (201 1 ) also identified a number of issues with existing 
satellite data such as the lack of a suitable atmospheric correction al- 
gorithm for multi-angular CHRIS/Proba data and a poor geo-location 
of the instrument. Our previous work has shown that some of these 
challenges could be overcome (Hilker et al., 2009b; Lyapustin et al., 
2007); and Hall et al. (2011) showed that despite these limitations 
a meaningful retrieval of s is still possible. However, at this current 
stage, the limitations existing with CHRIS/Proba data prevent it 
from being used routinely. Consequently, the approach introduced 
in this study needs to be seen more as demonstration and as a moti- 
vation for developing a MISR-like space mission to enable direct space 
borne observations of s; used in a data assimilation scheme to obtain 
spatially contiguous and temporally continuous estimates of vegeta- 
tion photosynthesis. 

While CHRIS/Proba data cannot be used in an operational sense, 
this satellite platform provides a unique opportunity to further test 
and develop the technique described in this study. At this point in 
its development we feel that the approach may be mature enough 
to study the relationship among vegetation light use efficiency, cli- 
mate and edaphic factors over transects with variations in these fac- 
tors. A new, multi-angular satellite sensor, as proposed in Hall et al. 
(2011), which uses 5 carefully selected wavebands to obtain f par 
from NDV1 and e from PRI (plus one band for atmospheric correction) 
would allow us to apply a technique such as applied in this study 
globally. In addition, if satellite data was available on a regular basis, 
this would also overcome the limitations of inferring the shape pa- 
rameters of the response curves from flux tower observations as 
they too could be directly obtained from space. 

5. Conclusion 

This study has demonstrated that spatially contiguous and tempo- 
rally continuous estimates of s can yield dramatically enhanced esti- 
mates of plant photosynthesis. We have shown that s varies greatly 
in both space and time and therefore requires comprehensive model- 
ing and sophisticated measurement techniques. The availability of 
multi-angular satellite data to infer instantaneous e on a routinely 
basis from space would allow us obtain more realistic estimates of 
vegetation carbon uptake through data assimilation and would there- 
fore be a significant contribution to closing the terrestrial carbon 
cycle. 
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